R Notebook

rm(list = ls())
library(sf)
link='https://github.com/MAGALLANESJoseManuel/Midis_Bono600/raw/main/midisceplanmap.geojson'
midis=sf::read_sf(link)

names(midis)

 [1] "DEPARTAMENTO"                  "PROVINCIA"                    
 [3] "DISTRITO"                      "INSTITUCIO"                   
 [5] "UBIGEO_yana"                   "UBIGEO"                       
 [7] "programado_yana"               "ejecutado_yana"               
 [9] "yanaProp"                      "PeruLibre_congre"             
[11] "congreEmitidos"                "proporcion_congre"            
[13] "PeruLibre_presi"               "presiEmitidos"                
[15] "proporcion_presi"              "programado_b600"              
[17] "ejecutado_b600"                "nivel"                        
[19] "pobTotal"                      "pobDiscapacidad"              
[21] "superficie"                    "densidad"                     
[23] "altitud"                       "latitud"                      
[25] "longitud"                      "numeroCentrosPoblados"        
[27] "porcentajeDesnutCronica5menos" "conteoPobres"                 
[29] "conteoPobresExtremos"          "porcentajePobrezaExtrema"     
[31] "b600prop"                      "geometry"                     

str(midis)

sf [1,874 × 32] (S3: sf/tbl_df/tbl/data.frame)
 $ DEPARTAMENTO                 : chr [1:1874] "TACNA" "TACNA" "TACNA" "TACNA" ...
 $ PROVINCIA                    : chr [1:1874] "TACNA" "TACNA" "TACNA" "TACNA" ...
 $ DISTRITO                     : chr [1:1874] "CORONEL GREGORIO ALBARRACIN LANCHIPA" "POCOLLAY" "CALANA" "TACNA" ...
 $ INSTITUCIO                   : chr [1:1874] "IGN" "IGN" "IGN" "IGN" ...
 $ UBIGEO_yana                  : int [1:1874] 230110 230108 230103 230101 230109 230104 230303 180301 230102 230106 ...
 $ UBIGEO                       : int [1:1874] 230110 230108 230103 230101 230109 230104 230303 180301 230102 230106 ...
 $ programado_yana              : int [1:1874] 50230 7204 1879 40771 1766 21048 1076 26281 17153 1227 ...
 $ ejecutado_yana               : int [1:1874] 46451 6746 1753 37323 1625 19352 973 24819 15743 1133 ...
 $ yanaProp                     : num [1:1874] 0.925 0.936 0.933 0.915 0.92 ...
 $ PeruLibre_congre             : num [1:1874] 12877 1776 622 8565 353 ...
 $ congreEmitidos               : num [1:1874] 58524 10919 3117 72366 1664 ...
 $ proporcion_congre            : num [1:1874] 0.22 0.163 0.2 0.118 0.212 ...
 $ PeruLibre_presi              : num [1:1874] 19737 2845 943 12927 547 ...
 $ presiEmitidos                : num [1:1874] 59715 11168 3117 73341 1664 ...
 $ proporcion_presi             : num [1:1874] 0.331 0.255 0.303 0.176 0.329 ...
 $ programado_b600              : int [1:1874] 28770 4445 1218 26312 1240 11915 0 15989 10685 720 ...
 $ ejecutado_b600               : int [1:1874] 26662 4130 1158 24021 1171 10882 0 14975 9794 681 ...
 $ nivel                        : chr [1:1874] "Distrito" "Distrito" "Distrito" "Distrito" ...
 $ pobTotal                     : num [1:1874] 104353 17856 4441 130965 2747 ...
 $ pobDiscapacidad              : num [1:1874] 953 165 31 1436 9 ...
 $ superficie                   : num [1:1874] 188 266 108 1878 1116 ...
 $ densidad                     : num [1:1874] 555.84 67.22 40.98 69.74 2.46 ...
 $ altitud                      : num [1:1874] 562 690 881 583 397 ...
 $ latitud                      : num [1:1874] -18 -18 -17.9 -18 -17.9 ...
 $ longitud                     : num [1:1874] -70.3 -70.2 -70.2 -70.3 -70.6 ...
 $ numeroCentrosPoblados        : num [1:1874] 32 13 37 18 33 15 17 4 8 22 ...
 $ porcentajeDesnutCronica5menos: num [1:1874] 3.3 2.81 3.33 2.76 3.73 ...
 $ conteoPobres                 : int [1:1874] 19552 2393 324 10072 292 8185 235 4806 5842 418 ...
 $ conteoPobresExtremos         : int [1:1874] 1715 168 19 296 31 536 16 241 368 42 ...
 $ porcentajePobrezaExtrema     : num [1:1874] 1.644 0.944 0.432 0.226 1.143 ...
 $ b600prop                     : num [1:1874] 0.927 0.929 0.951 0.913 0.944 ...
 $ geometry                     :sfc_MULTIPOLYGON of length 1874; first list element: List of 1
  ..$ :List of 1
  .. ..$ : num [1:282, 1:2] -70.2 -70.2 -70.2 -70.2 -70.2 ...
  ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
  ..- attr(*, "names")= chr [1:31] "DEPARTAMENTO" "PROVINCIA" "DISTRITO" "INSTITUCIO" ...

midis[,c('ejecutado_yana','conteoPobres')]

Simple feature collection with 1874 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.32816 ymin: -18.35067 xmax: -68.65218 ymax: -0.03877659
Geodetic CRS:  WGS 84

h1=formula(ejecutado_yana~proporcion_presi+proporcion_congre+porcentajePobrezaExtrema)
reg1=glm(h1,offset = log(pobTotal),family = 'poisson',data=na.omit(midis))
summary(reg1)

Call:
glm(formula = h1, family = "poisson", data = na.omit(midis), 
    offset = log(pobTotal))

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-195.792    -3.013     5.420    12.695    95.230  

Coefficients:
                           Estimate Std. Error  z value Pr(>|z|)    
(Intercept)              -1.354e+00  5.941e-04 -2278.86   <2e-16 ***
proporcion_presi          1.703e+00  1.175e-02   144.90   <2e-16 ***
proporcion_congre        -6.546e-01  1.739e-02   -37.65   <2e-16 ***
porcentajePobrezaExtrema  1.222e-02  4.753e-05   257.12   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 875774  on 887  degrees of freedom
Residual deviance: 445901  on 884  degrees of freedom
AIC: 454507

Number of Fisher Scoring iterations: 4

exp(coef(reg1))

             (Intercept)         proporcion_presi        proporcion_congre 
               0.2582633                5.4902408                0.5196249 
porcentajePobrezaExtrema 
               1.0122956 

library(ggplot2)
base=ggplot(data=midis,aes(x=b600prop,y=yanaProp)) 
base + geom_point()

sd(midis$ejecutado_yana)/mean(midis$ejecutado_yana)

midis$propDisca=midis$pobDiscapacidad/midis$pobTotal

AER::dispersiontest(reg1,alternative='greater')$ p.value<0.05

[1] TRUE

h2=formula(ejecutado_yana~proporcion_presi+proporcion_congre+ porcentajeDesnutCronica5menos+numeroCentrosPoblados+offset(log(programado_yana)))

#na.omit(midis)
rbn=MASS::glm.nb(h2,data=na.omit(midis))#,control=glm.control(maxit=50))

summary(rbn)

Call:
MASS::glm.nb(formula = h2, data = na.omit(midis), init.theta = 1158.09909, 
    link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-9.0667  -0.2518   0.1661   0.5200   1.7103  

Coefficients:
                                Estimate Std. Error z value Pr(>|z|)    
(Intercept)                   -7.841e-02  3.020e-03 -25.961  < 2e-16 ***
proporcion_presi               1.762e-01  2.808e-02   6.274 3.52e-10 ***
proporcion_congre             -1.452e-01  3.853e-02  -3.768 0.000165 ***
porcentajeDesnutCronica5menos -1.439e-04  1.434e-04  -1.004 0.315613    
numeroCentrosPoblados         -8.471e-05  2.310e-05  -3.667 0.000246 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(1158.099) family taken to be 1)

    Null deviance: 808.82  on 887  degrees of freedom
Residual deviance: 728.85  on 883  degrees of freedom
AIC: 10518

Number of Fisher Scoring iterations: 1

              Theta:  1158.1 
          Std. Err.:  71.7 

 2 x log-likelihood:  -10505.95 

rbn=MASS::glm.nb(h2,data=midis,control=glm.control(maxit=50))

summary(rbn)

exp(coef(rbn))

                  (Intercept)              proporcion_presi 
                    0.3688043                     1.4071911 
            proporcion_congre porcentajeDesnutCronica5menos 
                    1.2191827                     1.0071285 
        numeroCentrosPoblados 
                    1.0000267 

library(ggplot2)
dotwhisker::dwplot(rbn,exp=T)+ geom_vline(
           xintercept = 1,
           colour = "grey60",
           linetype = 2
       )
# + scale_y_discrete(labels=c("trabajo\nindependiente","analfabetismo")) + scale_color_discrete(name="Modelos para:\nCantidad de Asegurados")